home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / schur < prev    next >
Text File  |  1994-03-17  |  18KB  |  669 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*    
  28.     File containing routines for computing the Schur decomposition
  29.     of a real non-symmetric matrix
  30.     See also: hessen.c
  31. */
  32.  
  33. #include    <stdio.h>
  34. #include    <math.h>
  35. #include    "matrix.h"
  36. #include        "matrix2.h"
  37.  
  38.  
  39. static char rcsid[] = "$Id: schur.c,v 1.7 1994/03/17 05:36:53 des Exp $";
  40.  
  41.  
  42.  
  43. #ifndef ANSI_C
  44. static    void    hhldr3(x,y,z,nu1,beta,newval)
  45. double    x, y, z;
  46. Real    *nu1, *beta, *newval;
  47. #else
  48. static    void    hhldr3(double x, double y, double z,
  49.                Real *nu1, Real *beta, Real *newval)
  50. #endif
  51. {
  52.     Real    alpha;
  53.  
  54.     if ( x >= 0.0 )
  55.         alpha = sqrt(x*x+y*y+z*z);
  56.     else
  57.         alpha = -sqrt(x*x+y*y+z*z);
  58.     *nu1 = x + alpha;
  59.     *beta = 1.0/(alpha*(*nu1));
  60.     *newval = alpha;
  61. }
  62.  
  63. #ifndef ANSI_C
  64. static    void    hhldr3cols(A,k,j0,beta,nu1,nu2,nu3)
  65. MAT    *A;
  66. int    k, j0;
  67. double    beta, nu1, nu2, nu3;
  68. #else
  69. static    void    hhldr3cols(MAT *A, int k, int j0, double beta,
  70.                double nu1, double nu2, double nu3)
  71. #endif
  72. {
  73.     Real    **A_me, ip, prod;
  74.     int    j, n;
  75.  
  76.     if ( k < 0 || k+3 > A->m || j0 < 0 )
  77.         error(E_BOUNDS,"hhldr3cols");
  78.     A_me = A->me;        n = A->n;
  79.  
  80.     /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, A at 0x%lx, m = %d, n = %d\n",
  81.            __LINE__, j0, k, (long)A, A->m, A->n); */
  82.     /* printf("hhldr3cols: A (dumped) =\n");    m_dump(stdout,A); */
  83.  
  84.     for ( j = j0; j < n; j++ )
  85.     {
  86.         /*****        
  87.         ip = nu1*A_me[k][j] + nu2*A_me[k+1][j] + nu3*A_me[k+2][j];
  88.         prod = ip*beta;
  89.         A_me[k][j]   -= prod*nu1;
  90.         A_me[k+1][j] -= prod*nu2;
  91.         A_me[k+2][j] -= prod*nu3;
  92.         *****/
  93.         /* printf("hhldr3cols: j = %d\n", j); */
  94.  
  95.         ip = nu1*m_entry(A,k,j)+nu2*m_entry(A,k+1,j)+nu3*m_entry(A,k+2,j);
  96.         prod = ip*beta;
  97.         /*****
  98.         m_set_val(A,k  ,j,m_entry(A,k  ,j) - prod*nu1);
  99.         m_set_val(A,k+1,j,m_entry(A,k+1,j) - prod*nu2);
  100.         m_set_val(A,k+2,j,m_entry(A,k+2,j) - prod*nu3);
  101.         *****/
  102.         m_add_val(A,k  ,j,-prod*nu1);
  103.         m_add_val(A,k+1,j,-prod*nu2);
  104.         m_add_val(A,k+2,j,-prod*nu3);
  105.  
  106.     }
  107.     /* printf("hhldr3cols:(l.%d) j0 = %d, k = %d, m = %d, n = %d\n",
  108.            __LINE__, j0, k, A->m, A->n); */
  109.     /* putc('\n',stdout); */
  110. }
  111.  
  112. #ifndef ANSI_C
  113. static    void    hhldr3rows(A,k,i0,beta,nu1,nu2,nu3)
  114. MAT    *A;
  115. int    k, i0;
  116. double    beta, nu1, nu2, nu3;
  117. #else
  118. static    void    hhldr3rows(MAT *A, int k, int i0, double beta,
  119.                double nu1, double nu2, double nu3)
  120. #endif
  121. {
  122.     Real    **A_me, ip, prod;
  123.     int    i, m;
  124.  
  125.     /* printf("hhldr3rows:(l.%d) A at 0x%lx\n", __LINE__, (long)A); */
  126.     /* printf("hhldr3rows: k = %d\n", k); */
  127.     if ( k < 0 || k+3 > A->n )
  128.         error(E_BOUNDS,"hhldr3rows");
  129.     A_me = A->me;        m = A->m;
  130.     i0 = min(i0,m-1);
  131.  
  132.     for ( i = 0; i <= i0; i++ )
  133.     {
  134.         /****
  135.         ip = nu1*A_me[i][k] + nu2*A_me[i][k+1] + nu3*A_me[i][k+2];
  136.         prod = ip*beta;
  137.         A_me[i][k]   -= prod*nu1;
  138.         A_me[i][k+1] -= prod*nu2;
  139.         A_me[i][k+2] -= prod*nu3;
  140.         ****/
  141.  
  142.         ip = nu1*m_entry(A,i,k)+nu2*m_entry(A,i,k+1)+nu3*m_entry(A,i,k+2);
  143.         prod = ip*beta;
  144.         m_add_val(A,i,k  , - prod*nu1);
  145.         m_add_val(A,i,k+1, - prod*nu2);
  146.         m_add_val(A,i,k+2, - prod*nu3);
  147.  
  148.     }
  149. }
  150.  
  151. /* schur -- computes the Schur decomposition of the matrix A in situ
  152.     -- optionally, gives Q matrix such that Q^T.A.Q is upper triangular
  153.     -- returns upper triangular Schur matrix */
  154. MAT    *schur(A,Q)
  155. MAT    *A, *Q;
  156. {
  157.     int        i, j, iter, k, k_min, k_max, k_tmp, n, split;
  158.     Real    beta2, c, discrim, dummy, nu1, s, t, tmp, x, y, z;
  159.     Real    **A_me;
  160.     Real    sqrt_macheps;
  161.     static    VEC    *diag=VNULL, *beta=VNULL;
  162.     
  163.     if ( ! A )
  164.     error(E_NULL,"schur");
  165.     if ( A->m != A->n || ( Q && Q->m != Q->n ) )
  166.     error(E_SQUARE,"schur");
  167.     if ( Q != MNULL && Q->m != A->m )
  168.     error(E_SIZES,"schur");
  169.     n = A->n;
  170.     diag = v_resize(diag,A->n);
  171.     beta = v_resize(beta,A->n);
  172.     MEM_STAT_REG(diag,TYPE_VEC);
  173.     MEM_STAT_REG(beta,TYPE_VEC);
  174.     /* compute Hessenberg form */
  175.     Hfactor(A,diag,beta);
  176.     
  177.     /* save Q if necessary */
  178.     if ( Q )
  179.     Q = makeHQ(A,diag,beta,Q);
  180.     makeH(A,A);
  181.  
  182.     sqrt_macheps = sqrt(MACHEPS);
  183.  
  184.     k_min = 0;    A_me = A->me;
  185.  
  186.     while ( k_min < n )
  187.     {
  188.     Real    a00, a01, a10, a11;
  189.     double    scale, t, numer, denom;
  190.  
  191.     /* find k_max to suit:
  192.        submatrix k_min..k_max should be irreducible */
  193.     k_max = n-1;
  194.     for ( k = k_min; k < k_max; k++ )
  195.         /* if ( A_me[k+1][k] == 0.0 ) */
  196.         if ( m_entry(A,k+1,k) == 0.0 )
  197.         {    k_max = k;    break;    }
  198.  
  199.     if ( k_max <= k_min )
  200.     {
  201.         k_min = k_max + 1;
  202.         continue;        /* outer loop */
  203.     }
  204.  
  205.     /* check to see if we have a 2 x 2 block
  206.        with complex eigenvalues */
  207.     if ( k_max == k_min + 1 )
  208.     {
  209.         /* tmp = A_me[k_min][k_min] - A_me[k_max][k_max]; */
  210.         a00 = m_entry(A,k_min,k_min);
  211.         a01 = m_entry(A,k_min,k_max);
  212.         a10 = m_entry(A,k_max,k_min);
  213.         a11 = m_entry(A,k_max,k_max);
  214.         tmp = a00 - a11;
  215.         /* discrim = tmp*tmp +
  216.         4*A_me[k_min][k_max]*A_me[k_max][k_min]; */
  217.         discrim = tmp*tmp +
  218.         4*a01*a10;
  219.         if ( discrim < 0.0 )
  220.         {    /* yes -- e-vals are complex
  221.            -- put 2 x 2 block in form [a b; c a];
  222.            then eigenvalues have real part a & imag part sqrt(|bc|) */
  223.         numer = - tmp;
  224.         denom = ( a01+a10 >= 0.0 ) ?
  225.             (a01+a10) + sqrt((a01+a10)*(a01+a10)+tmp*tmp) :
  226.             (a01+a10) - sqrt((a01+a10)*(a01+a10)+tmp*tmp);
  227.         if ( denom != 0.0 )
  228.         {   /* t = s/c = numer/denom */
  229.             t = numer/denom;
  230.             scale = c = 1.0/sqrt(1+t*t);
  231.             s = c*t;
  232.         }
  233.         else
  234.         {
  235.             c = 1.0;
  236.             s = 0.0;
  237.         }
  238.         rot_cols(A,k_min,k_max,c,s,A);
  239.         rot_rows(A,k_min,k_max,c,s,A);
  240.         if ( Q != MNULL )
  241.             rot_cols(Q,k_min,k_max,c,s,Q);
  242.         k_min = k_max + 1;
  243.         continue;
  244.         }
  245.         else /* discrim >= 0; i.e. block has two real eigenvalues */
  246.         {    /* no -- e-vals are not complex;
  247.            split 2 x 2 block and continue */
  248.         /* s/c = numer/denom */
  249.         numer = ( tmp >= 0.0 ) ?
  250.             - tmp - sqrt(discrim) : - tmp + sqrt(discrim);
  251.         denom = 2*a01;
  252.         if ( fabs(numer) < fabs(denom) )
  253.         {   /* t = s/c = numer/denom */
  254.             t = numer/denom;
  255.             scale = c = 1.0/sqrt(1+t*t);
  256.             s = c*t;
  257.         }
  258.         else if ( numer != 0.0 )
  259.         {   /* t = c/s = denom/numer */
  260.             t = denom/numer;
  261.             scale = 1.0/sqrt(1+t*t);
  262.             c = fabs(t)*scale;
  263.             s = ( t >= 0.0 ) ? scale : -scale;
  264.         }
  265.         else /* numer == denom == 0 */
  266.         {
  267.             c = 0.0;
  268.             s = 1.0;
  269.         }
  270.         rot_cols(A,k_min,k_max,c,s,A);
  271.         rot_rows(A,k_min,k_max,c,s,A);
  272.         /* A->me[k_max][k_min] = 0.0; */
  273.         if ( Q != MNULL )
  274.             rot_cols(Q,k_min,k_max,c,s,Q);
  275.         k_min = k_max + 1;    /* go to next block */
  276.         continue;
  277.         }
  278.     }
  279.  
  280.     /* now have r x r block with r >= 2:
  281.        apply Francis QR step until block splits */
  282.     split = FALSE;        iter = 0;
  283.     while ( ! split )
  284.     {
  285.         iter++;
  286.         
  287.         /* set up Wilkinson/Francis complex shift */
  288.         k_tmp = k_max - 1;
  289.  
  290.         a00 = m_entry(A,k_tmp,k_tmp);
  291.         a01 = m_entry(A,k_tmp,k_max);
  292.         a10 = m_entry(A,k_max,k_tmp);
  293.         a11 = m_entry(A,k_max,k_max);
  294.  
  295.         /* treat degenerate cases differently
  296.            -- if there are still no splits after five iterations
  297.               and the bottom 2 x 2 looks degenerate, force it to
  298.           split */
  299.         if ( iter >= 5 &&
  300.          fabs(a00-a11) < sqrt_macheps*(fabs(a00)+fabs(a11)) &&
  301.          (fabs(a01) < sqrt_macheps*(fabs(a00)+fabs(a11)) ||
  302.           fabs(a10) < sqrt_macheps*(fabs(a00)+fabs(a11))) )
  303.         {
  304.           if ( fabs(a01) < sqrt_macheps*(fabs(a00)+fabs(a11)) )
  305.         m_set_val(A,k_tmp,k_max,0.0);
  306.           if ( fabs(a10) < sqrt_macheps*(fabs(a00)+fabs(a11)) )
  307.         {
  308.           m_set_val(A,k_max,k_tmp,0.0);
  309.           split = TRUE;
  310.           continue;
  311.         }
  312.         }
  313.  
  314.         s = a00 + a11;
  315.         t = a00*a11 - a01*a10;
  316.  
  317.         /* break loop if a 2 x 2 complex block */
  318.         if ( k_max == k_min + 1 && s*s < 4.0*t )
  319.         {
  320.         split = TRUE;
  321.         continue;
  322.         }
  323.  
  324.         /* perturb shift if convergence is slow */
  325.         if ( (iter % 10) == 0 )
  326.         {    s += iter*0.02;        t += iter*0.02;
  327.         }
  328.  
  329.         /* set up Householder transformations */
  330.         k_tmp = k_min + 1;
  331.         /********************
  332.         x = A_me[k_min][k_min]*A_me[k_min][k_min] +
  333.         A_me[k_min][k_tmp]*A_me[k_tmp][k_min] -
  334.             s*A_me[k_min][k_min] + t;
  335.         y = A_me[k_tmp][k_min]*
  336.         (A_me[k_min][k_min]+A_me[k_tmp][k_tmp]-s);
  337.         if ( k_min + 2 <= k_max )
  338.         z = A_me[k_tmp][k_min]*A_me[k_min+2][k_tmp];
  339.         else
  340.         z = 0.0;
  341.         ********************/
  342.  
  343.         a00 = m_entry(A,k_min,k_min);
  344.         a01 = m_entry(A,k_min,k_tmp);
  345.         a10 = m_entry(A,k_tmp,k_min);
  346.         a11 = m_entry(A,k_tmp,k_tmp);
  347.  
  348.         /********************
  349.         a00 = A->me[k_min][k_min];
  350.         a01 = A->me[k_min][k_tmp];
  351.         a10 = A->me[k_tmp][k_min];
  352.         a11 = A->me[k_tmp][k_tmp];
  353.         ********************/
  354.         x = a00*a00 + a01*a10 - s*a00 + t;
  355.         y = a10*(a00+a11-s);
  356.         if ( k_min + 2 <= k_max )
  357.         z = a10* /* m_entry(A,k_min+2,k_tmp) */ A->me[k_min+2][k_tmp];
  358.         else
  359.         z = 0.0;
  360.  
  361.         for ( k = k_min; k <= k_max-1; k++ )
  362.         {
  363.         if ( k < k_max - 1 )
  364.         {
  365.             hhldr3(x,y,z,&nu1,&beta2,&dummy);
  366.             tracecatch(hhldr3cols(A,k,max(k-1,0),  beta2,nu1,y,z),"schur");
  367.             tracecatch(hhldr3rows(A,k,min(n-1,k+3),beta2,nu1,y,z),"schur");
  368.             if ( Q != MNULL )
  369.             hhldr3rows(Q,k,n-1,beta2,nu1,y,z);
  370.         }
  371.         else
  372.         {
  373.             givens(x,y,&c,&s);
  374.             rot_cols(A,k,k+1,c,s,A);
  375.             rot_rows(A,k,k+1,c,s,A);
  376.             if ( Q )
  377.             rot_cols(Q,k,k+1,c,s,Q);
  378.         }
  379.         /* if ( k >= 2 )
  380.             m_set_val(A,k,k-2,0.0); */
  381.         /* x = A_me[k+1][k]; */
  382.         x = m_entry(A,k+1,k);
  383.         if ( k <= k_max - 2 )
  384.             /* y = A_me[k+2][k];*/
  385.             y = m_entry(A,k+2,k);
  386.         else
  387.             y = 0.0;
  388.         if ( k <= k_max - 3 )
  389.             /* z = A_me[k+3][k]; */
  390.             z = m_entry(A,k+3,k);
  391.         else
  392.             z = 0.0;
  393.         }
  394.         /* if ( k_min > 0 )
  395.         m_set_val(A,k_min,k_min-1,0.0);
  396.         if ( k_max < n - 1 )
  397.         m_set_val(A,k_max+1,k_max,0.0); */
  398.         for ( k = k_min; k <= k_max-2; k++ )
  399.         {
  400.         /* zero appropriate sub-diagonals */
  401.         m_set_val(A,k+2,k,0.0);
  402.         if ( k < k_max-2 )
  403.             m_set_val(A,k+3,k,0.0);
  404.         }
  405.  
  406.         /* test to see if matrix should split */
  407.         for ( k = k_min; k < k_max; k++ )
  408.         if ( fabs(A_me[k+1][k]) < MACHEPS*
  409.             (fabs(A_me[k][k])+fabs(A_me[k+1][k+1])) )
  410.         {    A_me[k+1][k] = 0.0;    split = TRUE;    }
  411.     }
  412.     }
  413.     
  414.     /* polish up A by zeroing strictly lower triangular elements
  415.        and small sub-diagonal elements */
  416.     for ( i = 0; i < A->m; i++ )
  417.     for ( j = 0; j < i-1; j++ )
  418.         A_me[i][j] = 0.0;
  419.     for ( i = 0; i < A->m - 1; i++ )
  420.     if ( fabs(A_me[i+1][i]) < MACHEPS*
  421.         (fabs(A_me[i][i])+fabs(A_me[i+1][i+1])) )
  422.         A_me[i+1][i] = 0.0;
  423.  
  424.     return A;
  425. }
  426.  
  427. /* schur_vals -- compute real & imaginary parts of eigenvalues
  428.     -- assumes T contains a block upper triangular matrix
  429.         as produced by schur()
  430.     -- real parts stored in real_pt, imaginary parts in imag_pt */
  431. void    schur_evals(T,real_pt,imag_pt)
  432. MAT    *T;
  433. VEC    *real_pt, *imag_pt;
  434. {
  435.     int    i, n;
  436.     Real    discrim, **T_me;
  437.     Real    diff, sum, tmp;
  438.  
  439.     if ( ! T || ! real_pt || ! imag_pt )
  440.         error(E_NULL,"schur_evals");
  441.     if ( T->m != T->n )
  442.         error(E_SQUARE,"schur_evals");
  443.     n = T->n;    T_me = T->me;
  444.     real_pt = v_resize(real_pt,(u_int)n);
  445.     imag_pt = v_resize(imag_pt,(u_int)n);
  446.  
  447.     i = 0;
  448.     while ( i < n )
  449.     {
  450.         if ( i < n-1 && T_me[i+1][i] != 0.0 )
  451.         {   /* should be a complex eigenvalue */
  452.             sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  453.             diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  454.             discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  455.             if ( discrim < 0.0 )
  456.             {    /* yes -- complex e-vals */
  457.             real_pt->ve[i] = real_pt->ve[i+1] = sum;
  458.             imag_pt->ve[i] = sqrt(-discrim);
  459.             imag_pt->ve[i+1] = - imag_pt->ve[i];
  460.             }
  461.             else
  462.             {    /* no -- actually both real */
  463.             tmp = sqrt(discrim);
  464.             real_pt->ve[i]   = sum + tmp;
  465.             real_pt->ve[i+1] = sum - tmp;
  466.             imag_pt->ve[i]   = imag_pt->ve[i+1] = 0.0;
  467.             }
  468.             i += 2;
  469.         }
  470.         else
  471.         {   /* real eigenvalue */
  472.             real_pt->ve[i] = T_me[i][i];
  473.             imag_pt->ve[i] = 0.0;
  474.             i++;
  475.         }
  476.     }
  477. }
  478.  
  479. /* schur_vecs -- returns eigenvectors computed from the real Schur
  480.         decomposition of a matrix
  481.     -- T is the block upper triangular Schur matrix
  482.     -- Q is the orthognal matrix where A = Q.T.Q^T
  483.     -- if Q is null, the eigenvectors of T are returned
  484.     -- X_re is the real part of the matrix of eigenvectors,
  485.         and X_im is the imaginary part of the matrix.
  486.     -- X_re is returned */
  487. MAT    *schur_vecs(T,Q,X_re,X_im)
  488. MAT    *T, *Q, *X_re, *X_im;
  489. {
  490.     int    i, j, limit;
  491.     Real    t11_re, t11_im, t12, t21, t22_re, t22_im;
  492.     Real    l_re, l_im, det_re, det_im, invdet_re, invdet_im,
  493.         val1_re, val1_im, val2_re, val2_im,
  494.         tmp_val1_re, tmp_val1_im, tmp_val2_re, tmp_val2_im, **T_me;
  495.     Real    sum, diff, discrim, magdet, norm, scale;
  496.     static VEC    *tmp1_re=VNULL, *tmp1_im=VNULL,
  497.             *tmp2_re=VNULL, *tmp2_im=VNULL;
  498.  
  499.     if ( ! T || ! X_re )
  500.         error(E_NULL,"schur_vecs");
  501.     if ( T->m != T->n || X_re->m != X_re->n ||
  502.         ( Q != MNULL && Q->m != Q->n ) ||
  503.         ( X_im != MNULL && X_im->m != X_im->n ) )
  504.         error(E_SQUARE,"schur_vecs");
  505.     if ( T->m != X_re->m ||
  506.         ( Q != MNULL && T->m != Q->m ) ||
  507.         ( X_im != MNULL && T->m != X_im->m ) )
  508.         error(E_SIZES,"schur_vecs");
  509.  
  510.     tmp1_re = v_resize(tmp1_re,T->m);
  511.     tmp1_im = v_resize(tmp1_im,T->m);
  512.     tmp2_re = v_resize(tmp2_re,T->m);
  513.     tmp2_im = v_resize(tmp2_im,T->m);
  514.     MEM_STAT_REG(tmp1_re,TYPE_VEC);
  515.     MEM_STAT_REG(tmp1_im,TYPE_VEC);
  516.     MEM_STAT_REG(tmp2_re,TYPE_VEC);
  517.     MEM_STAT_REG(tmp2_im,TYPE_VEC);
  518.  
  519.     T_me = T->me;
  520.     i = 0;
  521.     while ( i < T->m )
  522.     {
  523.         if ( i+1 < T->m && T->me[i+1][i] != 0.0 )
  524.         {    /* complex eigenvalue */
  525.         sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  526.         diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  527.         discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  528.         l_re = l_im = 0.0;
  529.         if ( discrim < 0.0 )
  530.         {    /* yes -- complex e-vals */
  531.             l_re = sum;
  532.             l_im = sqrt(-discrim);
  533.         }
  534.         else /* not correct Real Schur form */
  535.             error(E_RANGE,"schur_vecs");
  536.         }
  537.         else
  538.         {
  539.         l_re = T_me[i][i];
  540.         l_im = 0.0;
  541.         }
  542.  
  543.         v_zero(tmp1_im);
  544.         v_rand(tmp1_re);
  545.         sv_mlt(MACHEPS,tmp1_re,tmp1_re);
  546.  
  547.         /* solve (T-l.I)x = tmp1 */
  548.         limit = ( l_im != 0.0 ) ? i+1 : i;
  549.         /* printf("limit = %d\n",limit); */
  550.         for ( j = limit+1; j < T->m; j++ )
  551.         tmp1_re->ve[j] = 0.0;
  552.         j = limit;
  553.         while ( j >= 0 )
  554.         {
  555.         if ( j > 0 && T->me[j][j-1] != 0.0 )
  556.         {   /* 2 x 2 diagonal block */
  557.             /* printf("checkpoint A\n"); */
  558.             val1_re = tmp1_re->ve[j-1] -
  559.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  560.             /* printf("checkpoint B\n"); */
  561.             val1_im = tmp1_im->ve[j-1] -
  562.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  563.             /* printf("checkpoint C\n"); */
  564.             val2_re = tmp1_re->ve[j] -
  565.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  566.             /* printf("checkpoint D\n"); */
  567.             val2_im = tmp1_im->ve[j] -
  568.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  569.             /* printf("checkpoint E\n"); */
  570.             
  571.             t11_re = T_me[j-1][j-1] - l_re;
  572.             t11_im = - l_im;
  573.             t22_re = T_me[j][j] - l_re;
  574.             t22_im = - l_im;
  575.             t12 = T_me[j-1][j];
  576.             t21 = T_me[j][j-1];
  577.  
  578.             scale =  fabs(T_me[j-1][j-1]) + fabs(T_me[j][j]) +
  579.             fabs(t12) + fabs(t21) + fabs(l_re) + fabs(l_im);
  580.  
  581.             det_re = t11_re*t22_re - t11_im*t22_im - t12*t21;
  582.             det_im = t11_re*t22_im + t11_im*t22_re;
  583.             magdet = det_re*det_re+det_im*det_im;
  584.             if ( sqrt(magdet) < MACHEPS*scale )
  585.             {
  586.                 det_re = MACHEPS*scale;
  587.             magdet = det_re*det_re+det_im*det_im;
  588.             }
  589.             invdet_re =   det_re/magdet;
  590.             invdet_im = - det_im/magdet;
  591.             tmp_val1_re = t22_re*val1_re-t22_im*val1_im-t12*val2_re;
  592.             tmp_val1_im = t22_im*val1_re+t22_re*val1_im-t12*val2_im;
  593.             tmp_val2_re = t11_re*val2_re-t11_im*val2_im-t21*val1_re;
  594.             tmp_val2_im = t11_im*val2_re+t11_re*val2_im-t21*val1_im;
  595.             tmp1_re->ve[j-1] = invdet_re*tmp_val1_re -
  596.                     invdet_im*tmp_val1_im;
  597.             tmp1_im->ve[j-1] = invdet_im*tmp_val1_re +
  598.                     invdet_re*tmp_val1_im;
  599.             tmp1_re->ve[j]   = invdet_re*tmp_val2_re -
  600.                     invdet_im*tmp_val2_im;
  601.             tmp1_im->ve[j]   = invdet_im*tmp_val2_re +
  602.                     invdet_re*tmp_val2_im;
  603.             j -= 2;
  604.             }
  605.             else
  606.         {
  607.             t11_re = T_me[j][j] - l_re;
  608.             t11_im = - l_im;
  609.             magdet = t11_re*t11_re + t11_im*t11_im;
  610.             scale = fabs(T_me[j][j]) + fabs(l_re);
  611.             if ( sqrt(magdet) < MACHEPS*scale )
  612.             {
  613.                 t11_re = MACHEPS*scale;
  614.             magdet = t11_re*t11_re + t11_im*t11_im;
  615.             }
  616.             invdet_re =   t11_re/magdet;
  617.             invdet_im = - t11_im/magdet;
  618.             /* printf("checkpoint F\n"); */
  619.             val1_re = tmp1_re->ve[j] -
  620.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  621.             /* printf("checkpoint G\n"); */
  622.             val1_im = tmp1_im->ve[j] -
  623.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  624.             /* printf("checkpoint H\n"); */
  625.             tmp1_re->ve[j] = invdet_re*val1_re - invdet_im*val1_im;
  626.             tmp1_im->ve[j] = invdet_im*val1_re + invdet_re*val1_im;
  627.             j -= 1;
  628.         }
  629.         }
  630.  
  631.         norm = v_norm_inf(tmp1_re) + v_norm_inf(tmp1_im);
  632.         sv_mlt(1/norm,tmp1_re,tmp1_re);
  633.         if ( l_im != 0.0 )
  634.         sv_mlt(1/norm,tmp1_im,tmp1_im);
  635.         mv_mlt(Q,tmp1_re,tmp2_re);
  636.         if ( l_im != 0.0 )
  637.         mv_mlt(Q,tmp1_im,tmp2_im);
  638.         if ( l_im != 0.0 )
  639.         norm = sqrt(in_prod(tmp2_re,tmp2_re)+in_prod(tmp2_im,tmp2_im));
  640.         else
  641.         norm = v_norm2(tmp2_re);
  642.         sv_mlt(1/norm,tmp2_re,tmp2_re);
  643.         if ( l_im != 0.0 )
  644.         sv_mlt(1/norm,tmp2_im,tmp2_im);
  645.  
  646.         if ( l_im != 0.0 )
  647.         {
  648.         if ( ! X_im )
  649.         error(E_NULL,"schur_vecs");
  650.         set_col(X_re,i,tmp2_re);
  651.         set_col(X_im,i,tmp2_im);
  652.         sv_mlt(-1.0,tmp2_im,tmp2_im);
  653.         set_col(X_re,i+1,tmp2_re);
  654.         set_col(X_im,i+1,tmp2_im);
  655.         i += 2;
  656.         }
  657.         else
  658.         {
  659.         set_col(X_re,i,tmp2_re);
  660.         if ( X_im != MNULL )
  661.             set_col(X_im,i,tmp1_im);    /* zero vector */
  662.         i += 1;
  663.         }
  664.     }
  665.  
  666.     return X_re;
  667. }
  668.  
  669.